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ill 


Rlood routing, forming one of the important facets 
of flood studies, has achieved vast progress in the last 
two decades due to general availahility of digital computers. 
Ainien*s implicit method is one popular method in this 
category, Ihe finite Element Method (EEM) which is gaining 
popularity in various "branches of mechanics has a great 
possibility of producing an effective numerical computation 
technique for handling a variety of flood routing problems. 

In this thesis an algorithm for flood routing in natural 
channels through the use of EEM with an implicit solution 
procedure is presented. 

In this algorithm, one dimensional method of analysis 
is adopted and provision has been made for inclusion of 
lateral flov/. The banded property of the global matrix 
has been used in finding a fast solution. The performance 
of the EE&I algorithm is compared with that of itoien’ s 
implicit solution procedure. The computer GPU time, errors 
and stability aspects are studied. The PEM algorithm is 
found to give stable and convergent results at a very wide 
range of time steps. 
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OHi ^ ER - I 
IITROEUGTI OIT ^ 

The importance of floods in human activities has 
generated a large number of studies towards understanding 
the various aspects of the flood flow* Eloods are caused 
mainly by heavy rainfall in the catchment melting of snow 
or heavy releases at an upstream storage station. As the 
conveying channel may not be adequate to safely carry the 
incoming high discharges, some precautionary and other mea- 
-sures have to be taken in advance to prevent the Ion of 
life and property. It requires the advance knowledge 
of the flood volume, its peak time, maximum stage to 
which it is likely to go, etc. at the selected points. 

This information may also be helpful for the regulation 
of different hydraulic structures in the system. This 
process by which we find out the channel outflows at 
certain selected stations of the channel, by knowing 
the upstream inflov/s in the channel is called flood 
routing. Inflows and outflows are in the form of stage 
or discharge hydrographs. 

Stream flow is a phase of hydrologic cycle and is 
not well defined process. Actually it is a unification 
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of various physical processes. There is a continuous 
change from one state to another state with respect to 
space as well as time. For example parameters like 
channel geometry, hed roughness, stage etc. can get 
changed in space and time. Sventhough these can be 
modelled to a great extent, some assumptions and simpli- 
fications are necessary. One of the major simplification 
is that continuous processes are treated as discrete 
ones, v/holly or to some extent. In these, we should have 
a check that our models and its various assumptions 
and simplifications do not have adverse effects on the 
resulting output. 

Flood routing may be done by either a process 
approch or a system approach. If flood routing is done 
by using the process approach it is called hydraulic 
flood routing and by system approach, it is called 
hydrologic flood routing. The first one involve an 
mathematical modelling and the later one the conceptual 
modelling. These distinctions may be appropriate for 
well defined channels. In the case of natural channels, 
it may not be possible to have exact mathematical formu- 
lations and we have to apply conceptual approach at some 
places, hue to the remarkable development in the fields 



5 


of high speed digital computers and numerical methods, 
hydraulic models axe playing a significantly large role 
in this field. 

Unsteady flow in a rigid bed, open channel 
can be described by the St, Tenant equations. These 
equations are the equations of conservation of massand 
momentum and were first given by Saint Tenant in 
ISy-j ( > 5 ,8, 1 5) . Tenant equations are two non*- 

linear, first order, first degree partial differential 
equations. There are no mathematical solution to them 
and they can be solved numerically for specified boundary 
and Initial conditions. 

The solution procedures of St. Tenant equations 
are discusses in Chapter II. Chapter III discuss briefly 
the commonly used finite difference method (PUM) and its 
limitations, in EBM algorithm is developed and its 
characteristics enumerated in Chapter IT. A comparative 
study of the BUM and the BEM algorithm is presented 
with the help of two examples in Chapter T. The 
significant conclusions of this study are collated in 
Chapter TI, ippendix - A and B presents the listing 
of the BUM and BSM algorithms implemented through B0RTRAN'-10 
in the l)EC’-1090 system at Indian Institute of Technology, 
Kanpur. 



CHiPTSR - II 


RUIvIERI Gil SOLUTIO N OF T HE SiINI > -VSRiNT EQUi ll OJ^S_ 


2 • B ASIG EQTJjglORS 

The St. Yenant equations for unsteady flow are 

( 1 ) 

given a.s^ ^ ; 



+ A 


8Y 

Sx 
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( 2 . 1 ) 


and 


3 1 


+ Y 


6V 


8x 


+ g 


3x 


g( Sq** 


Sf) + 

A 


= 0 

(2.2) 


where Y 
A 
X 

y 

t 

s. 


= average velocity 
= cross-sectional flow area 
= distance in the flow direction 
= water depth 
= time 
= hed slope 

= friction slope 


q = lateral incoming discharge per unit length 

u = average component of inflow velocity in x 
direction 

g = gravitational acceleration. 
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ihe first equation is the continuity equation 


and the second one is the momentum equation, 
of these equations assumes that^ j 

Derivation 

(i) 

the flo\? is one dimensional 


( ii) 

the fluid is homogeneous and incompressible 

(iii) 

the channel bed is fixed 


(iv) 

the channel alignment is approximated 
straight line 

to be 

(v) 

the bottom slope is small 


(vi) 

the flow is gradually varied 


(vii) 

the water surface is horizontal across 
section 

a cross 

(viii) 

tliG volocity is consta^nt across a cross section 

( ix) 

the wind resistance is neglected. 



S OLUTIOIT 0? ST.Y EM'INT SQU ATIOFS 

The Yarious methods of solving St. Tenant 
equations can be broadly classified as (a) approximate 
numerical method, (b) complete numerical method. The 
approximate methods are based on drastically curtailed 
equation of motion. The complete numerical method obtains 
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the solution by St. Yenant eq.aations. The Direct Method, 
Pinite Element Method and Method of Oharacteri sties 
(Rectangular grid and Oharac ter i sties nodes) belong 
to this category. These are further classified in the 
implicit and explicit method. Subraraanya^ ^ has 
given general description of direct method and of method 


of characteristics. Different finite differencing 
schemes are available for these tv/o methods^"^ ’ 


Einite Element Method (EEIi'I) is of recent origin and 

( fi) 

is in developing stage* Rrobably Cooley and Moin'* ^ 
are the first one to propose an algorithm for EEM to 


solve St. Tenant equations. 


2 ^cqijplTOTS 

To start the computation, the initial values 
of unknowns are required. These are discharges and 
water depths at all the nodal points along the river. 
Generally one assumes that the flo?/ to be steady before 
the start of flood wave. Hence initial values may be 
computed by the solution of gradually varied steady flow 
equations, 

BOGHDiRY C OHDITIOIfS 

Boundary conditions are must for a solution 

( 10 ) 


of a system in space and time, Liggett and Woolhiser 
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specified the-t one condition at the upstream is needed 
if the flow is subcritical and t?/o conditions are needed 
at the upstream if the floy/ is supercritical. In addition 
for subcritical flow one condition at the down stream 
boundary also is needed* These may be in the form of 
specified depth or velocity at a time. If the flovj 
changes from subcritical to supercritical and vice versa, 
the boundary specifications may get changed. 

The upstream boundary condition is generally a 
discharge or stage hydro graph. The down stream boundary 
condition is given by a rating curve. It carries the 
assumption that this section is not subject to back- 
water effects from dovrastream regulation. Hov^ever if 
some flo’w restriction exist the computation should 
either be carried out with the downstream boundary 
located at this restriction, or a rating curve obtained 
from back water computations may be used. 

2*5 DYh^iiaia Mm KimiiTIO MODELS 

Dynamic waves propogate in two systems of chara- 
cteristics; in the upstream and downstream direction. 

St. Tenant equations describe the dynamic waves. Upstream 
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travel speed is given ly 


- V - tgJ (2.3) 

and downstream travel speed is given "by 

0^ = V + f~gy' (2.4) 

Kinematic waves posses only one system of 
characteristics. Their theoretical speed can be given 
by Kleitz-Seddon law*'^^^. 


G = 


dA 


(2.5) 


In a natural flood wave both types of wave movement 
are present. The bulk of flood wave moves substantially 
as a kinematic wave, while the dynamic wave fronts move 
in front and behind the main body of the flood wave. 
Observations supports the theory of kinematic waves, but 
the flood wave does not steepen as much as predicted 

(g) 

by it. Henderson^ ' showed that the magnitude of the 
deviations from kinematic wave behaviour depends 
mainly on the ratio between bed slope and wave slope 
and on Kroude Number. 

Kor dynamic wave the waiting curve is a looped 


one and i s given by 
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"Sy * J”' jT ^ 2 

az s^g ax s^g at 

where is the normal discharge and can he given by 
Ohezys, Blannings or any other emperical resistance 
equations^ . • 

For kinematic wave the discharge is equal to 
normal discharge and the rating curve, is a linear curve 
is given by 



Q Q]^ ••• (2,7) 

The models based on kinematic waves are called 
Kinematic or Approximate models, while the models 
based on Dynamic waves are called Dynamic or Complete 
models^ . 


2.6 CSHiiranii CHARAGTERISTIOS 


Channel geometry may be described in various 
forms, Weinmann^ has distinguished four broad 
ways. 

(i) Replacement of actual river by a uniform channel 
for a total length of reach, 

(ii) Replacement of actual river by a series of 
prismatic channels of different width. 
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(iii) Direct use of polygonal section. 

(iv) Stochastic generation of cross-section. 

Resistance properites .ore usually based on Chezy's 
or Manning’s formula. 'Ihe relevant resistance coeffi- 
cients may be taken either constant or as a function of 
position and depth. 

2 . 7 R0UTI R& 

These consist of the time and space increments 
and the different weighting coefficients. Yariation of 
these parameters may affect the efficiency and accuracy 
of the numerical procedure adopted and hence the proper 
selection is of these great importance. 



OHA PTER^ - III 

3 . 1 IDTTRODUG ijqi^ 

imong the direct method and method of characteri- 
stics, the finite difference method using an Implicit 
numerical procedure has been generally accepted as the 
most convenient method, imong the various implicit 
finite differencing schemes, imeiiis Implicit Scheme^ 
is considered to be the most efficient one^^^^, in 
this Chapter an Implicit EPM ilgorithm which is a 
modification of the imien’s work^”^’^^ is developed. 

This IS used as a basis for comparison of the efficiency 
of EEM algorithm presented in the subsequent chapter 
of this thesis, 

1>STA ILEI) RSSGRIRTIOR 

The basic equations are taken in the form used 
by imein and Rang in 1970^^^ as 


Ll + i 


Y ^ 

^ 1 

3t B 

8x 

9x 

B 


(3.1) 


6Y 




dj 

5x 



(3.2) 


and 
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where B 


dy 


Here u has heen 


assumed zero. 


The equations (3.1) and (3.2) are replaced by the 
algebric finite difference equations and then the solution 
of these equations will be obtained. To find the numerical 
solution of equations (3.1) and (3.2) by implicit method 
the (x, t) plana is discretized in a rectangular net 
covering the whole river reach and the time for which 
computations are to be carried. The x axis will be 



l'IG.3.1 (x,t) PLiiNB IN MmilVS IMPLICIT SCHEME 

be representing the initial condition and the line 

paralled to it ?;ill denote times. (Pig. 3.1). The 

spacing between these lines will be ^ x. Similarly 


t is the time increment and t axis denotes the 



13 


upstream boundary. The Line, parallel to it, drawn 
at X = L, where 1 is the length of the reach will 
denote downstream boiindary. Line in between will 
represent different nodal stations (fig. 3.1). The 
forming grid will be double subscripted one. The values 
of X and.-rOi t need not to be constant. But however 
these has been taken constant in the present study. 

Assuming that all the variables are known at all 

nodes on the row having t = t^ (I’Q) and variables at 

inh i 

the row having t = t (RS) are unknovm where 
t^"^^ = t^ + ill:, t . In a four point method the values 
of any function and its derivatives with respect to time 
and space are expressed in terms of the values of 
function at the four grid points. If the grid is 
found by the state lines t = t, t=t^'^\ x=x^ 
and X = x^_^^ (PQ, RS, PR and QS in Pig. 3.1) and it is 
assumed that a point M (Pig. 3.1) is centered between 
the space grid line and positioned at t = t^ 

+ e(t^''’^ *- t^) where 0 is a weighting factor, then 
the function ( ) and its derivatives will be 

eCkf + + (1 - e)(ki+ kj^p 

2 


(3.3) 



u 


Sk 

9x 

and 

ilS 

6t 



(3.4) 


(3.5) 


Tha value of 6 may be assumed to be between 0.5 

and 1.0. Por a value of 1 the scheme is known as the 

fully implicit scheme and for a value of 0.5, it is 

known as the box scheme, imain and Pang in 1970*-"’^ have 

used box scheme and imein and Ohu in 1975 ^^^ have used 

the fully implicit scheme. As per Quinn and Wylie^”*'^^ 

implicit schemes are weakly stable for 6 in between p0,5 

a.nd 0.6. Chaudri and contractor have shown that the 

implicit values are more accurate but generally unsta,ble 

for a value close to 0.5. Converse is the case for a 

( 18 ") 

value of 0 = 1,0, woinmann^ took the value of 0 = 1. 
However in the present study various values of 9 between 
0,5 and 1,0 have been tested to knov/. the sensitivity 
of the scheme and to suggest the optimum value of 6 for 
use. 

By putting values of functions and variables i. e. 

y, V, , h, etc. obtained from equations 3.3, 3,4 

S t Qx 
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and. 3.5 in equations 3.1 and 3.2, finite difference 
form of these equations are obtained. In this friction 
slope at a point (i, j) is computed hy pfenning's 

formula 





(n^ ) 




(rJ 



(3.6) 


where rJ = / ?j (hydraulic radius) 

and = wetted parameter 

Hereafter time dependent variable will not carry 
any superscript, if time line is j+l . 


Pinite Difference form are given as 


Jj. (yi. q. a + a' + “V 

4A X 

I 

((c + c')(b+b') + (h+h') +(e+e'))+0.5 tq 
(w + w') = 0 (3.7) 


and 


%^^i’ "^i’ ^i+1’ "'^i+l^ 

= (a' ’ + a' ' ' + (h' '+b' ’ *)+ T ■ (c ' ’+c ' » ') 

gz."'. t 

(d’'+d' '')+ -- -(h' '+h* ' ')+ q(o’'+o'*') = 0 

2 2g 


(3.8) 
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where 


= ^1+1 + ?! 


b = 

2e(Y^ + Y 

i+1 ^ 

c = 

20 

(^i+i 


e = 

20 

^ ^i+ 1 


h = 

20 

(1. 

h^i 



% 

%+i 


■) 


w = 20 ( 


% 


1 

a+1 


+ — . -) 
B. 


= -7? , 
■^1+1 


G ' 


= 20 (yJ + 

= 20 (7^1 - yi ) 


e' = 20 - 

aI 

h' = 20 ( ~ + 

3? 

1 


) 


4 


+1 


®i+1 


( 3 . 9 ) 


? T 


h- = + Y.^1 


c ” = b 
d * ' = e 
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H’ • = 20( Sfi + Sf 


o ' ' =20 ( -=ii- + ™JiL_ ) 




■^+1 

a' ' ’ = 

c ^ 


b' • ' = 

-Y- «- Y^ . 
X. 1+1 


c ' » ' = 

b' 


d1 • » = 

e ’ 


h' ' ’ = 

20(s|j_ + 

i+1^ " ^ 

0 ' ' ' = 

20 ( 

i 

Y? , 

+ ^1.11 
■^+1 


where 0 = 1-9 


It amy he seen that the constants given in Eq. 3.9 
carrying a(’")are independent of the variahles at time 
(j+ 1 ). Similarly unscripted constants and the constants 
carrying (*’) are independent of the variable at time j, 
This property should be kept in mind to avoid duplicate 
computations of these constants. 

In equations 3 . 7 , 3 . 8 , and 3.9 only the variables 
with the superscript 3+I are the unloiowns. These are 
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^ (fT- 1 ) interior points provide 

equations for the 2- unknowns. . 3.^0 additional 
are available from the boundary conditions. 


equations 


If the stage hydograph at up stream boundary is 
given then 


&o (yi) = 71 - fi(t^+'’) = 0 
and If discharge hydro graph is given then 

&o ^1) = ^1 Ai - f^(t 3 +'’) = 0 (3,11; 

^here fi(t'' ) is stage or discharge value taken from 
the respective hydrograph as tte case may be. 

Similarly one equation will be obtained from down 
stream bound, ary condition. This equation will be 

% - Q (3.12) 

where Q will be calculated from equations 2.6 or 2.7 for 
dynamic and kinematic cases respectively, in the present 
study Normal flow, Q^, is calculated by using Manning’s 
equations. 

Now our equations constitute a set of 2N nonlinear 
algebric equations in 2N unknowns. The equations can be 
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assembled as follo? 7 S 

(y-,, Ti) =0 

^21 = JiCyi, Ti, y^^,, y.^^) 0 

for i = i to lf-1 (3.13) 

^(2i+1) = y.^^, y.^^) = q 

for i = 1 to iT -1 

^2n = % (y^-, Yjj) = 0 


Altbough sq. 3,13 carries 2K unknowns yet each 
equation contains a maximum of four unknowns and the 
whole matrix can come in banded matrix of width five. The 
generalized Ifewton iteration method is used herein to 
find the solution of this system. Relating residuals of 

Eq, 3.13 to the partial derivatives, as per lewton 
iteration method^ 


6 

6y. 


o 


9G„ 


^1 




dY. 




3Rj_ 

8yi 


. + 

yi 


dp. 


dy- 


dY. 


8p. 


"^i+l * 


9Pj_ 

3Y- ^ 
1+1 


^^i+1 


^2i 


for i = 1 to R-l 
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1!e 



4- 






” ^2IT 


(3.H) 


''‘Vhicli constitute a set of 2IT linear equations in 
2]! unknowns and in wkicb and dY^ for i having 

value from 1 to iT are the additive ■ corrections to he 
applied to the respective y^ and Yj_ values at the (,j+1) 
time. Knov/ing the values of dy^ and dY^ from the 
solution of above given set of equations the values of 
yj_ and Y^^ may be modified in this manner. The method 
used here is Gaussian elemination method using banded 
matrix technique. At the terminal iteration, when the 
values of dy^ and/or dY^ are less then the permissible 
errors, we take the values of variables found as final 
and advance for time step (,j+2). 

The expression for the partial derivatives are 
obtained by differentiating Sq. 3.7 and 3.8 partially 
with respect to the variables 


= 1 - ~ ((b + b') + (( e + e') 

Sy^i^ 2, iv X 


0(1 


Sf 




®i 


.(p).e 


(3.15) 
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2I‘. 


e(i 


^®i+1 ^ 




1-4.. 1 (^). , 0 

(B. .)2 ^ 1+1 

^1+1 ^ 


6F^ 


— — (c + c' ■- b. >- h.’) 

X 


(3.16) 

(3.17) 


8 5’-; ^ 4- 

1- — , 2s. ij 


67 


i+1 


4 X 


(c + c ' + h + h') 


(3.18) 


SG. 


QJh 


2e (-1 . 1 ^ X s^. ( ), 


h 


3. 

4f- > 


q X 

’ 2g 


7. 

1 


B 


) 

2 


) 


(3.19) 


6G. 

3y 


i+1 


29 (1 + £ A X S.. ^ 

3 


^ ^fi+1 (l5r~ (fj )i+i 


^4^ X 

- q. -=”— 7,- 

2g 


B, 


i+1 


i+1 


(%+l) 


®i+1 


) 


h+1 

(3.20) 


8_G^ 

eV. 


g ^ ^ ^ ~^h " - c” ' >- 49 7^^) 


^fi 7. 

Vi V" 


+ 20 X 


(3.21) 
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Q G. 

_ 1 

Q Y- . 
1+1 


^ /S X , 0 , 2 X 

g A t 2g 1+1^ 


+ 26 ^cxX 


S4? 


ill + Y. . (3.22) 

g 


^'^i+l 


Similiarl.};" partial derivatives has been found for 
boundaries. if the stage hydrograph is given at upstream 
the partial derivatives will be 



(5.23) 


9G, 


ay. 


= 0 


(5.24) 


and if the discharge hydro graph is given at the upstream 
the partial derivatives will be 


^ - V-1 (3.25) 

9y>i 

8 Gq 

~~ = - A. (3.26) 

8 Yi 

Similarly for Kinematic wave condition at d^vm 
stream boundary the values will be given by 
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6F 


IT 




= \ - I 




^Y^ 


IT 


~ %T 


2/3 


■ep ^5/3 


“ 8y IJ - 2-5 %)S^ 

(3.27) 

(3.28) 


and for dynamic wave condition the down strPam boundary 
these values will he given, by 


= Y B + -1. 2 2/3 3P s5/3 

3y.j ^ ^ %T 3 W ^IT 

- 2.5 sy^ ^ 0.5 




(3.29) 


• So,A ^ 


and 


8Y- 


At ^ + 0.5 


'^IT 


IT 




.(2V,, - Y. 


A X 


IT ''11^1' 


) 


A t 

(3.30) 


where 


\ ■ tt; 




%( \t - ^ 

1 2/3 ^/2 

H” % ^ 


t 


O SA X 

^^^V2 

(3.31) 

(3.32) 


- V>) 


and nj^ = Mannings roughness constant 
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3 . 3 HlQSRjjffinNG 

The above implicit method of numerical solution was 
programmed in FORTRilT-IO and implemented at DEO-IOBO at 
Indian Institute of Technology, Kanpur, Initial condition, 

of boundary condition, Po. of nodes on x and t axis, 
time increament, space increment, lateral flow, bed slope 
and the permissible termination errors are being read as the 
input data. However input has also to be given for the 
form of output needed viz. the stations v/here the output 
is needed and the fime increament after which out put is 
needed- Provision has been made for giving stage and 
discharge hydrograph. If the programme does not iterate 
to a desired level in a specified num.ber of iterations, 
it is mentioned in the output and the programme is 
termin ated. The variable viz. Area, Breadth, Parameter, 
Manning’s n who are time and space dependent are to be 
defined in the form of subrouting functions. There 
derivatives are also introduced in the same way. Upstream 
hydrograph, either stage on discharge, is given in a 
functional subroutine as a function of time. Lateral 
flow is available by subroutine which gives the 
discharge and their velocity component in the direction 
of flow at a particular time for all the lateral flows. 
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Eqs, 3.13 and from 3.15 to 3.22 has been given in the 
form of functions. All the details has been provided in 
the main body of the programme. Flo?; chart of the progr-amme 
is shown in Fig. 3.2. a listing of the FORTRAN Program 
(FDM) is included in Appendix - A. 

^ ^ CHARACTERISTlas OF FDM PRO GRill 

The size of distance steps does not have much 
influence on the solution accuracy, as long as the 
variation of the flow characterises along the channel 
can be appropriately represented! In the present study 
selection of distance steps has been done by the 
consideration of stability. In the method of characteristics 
selection of ^t is done by follov/ing the Courand criteria 
A i; 

( |y + Q j < -i). However in the present study even 

bv X “ 

greater ( oc 2.5 times Courant value) values produced an 
stable result. The 6 should be selected for giving 
accurate, stable and oonvervent results. In the present 
study it has been found that 6 = 0.6 satisfies the above 
properties to a great extent and consume less CPU time. 
However 0 = 0.55 - 0.65 gives acceptable results and it 
may be selected, in a way that it reproduces the desired 
hydro graph. 



26 


w 


I read ROUTim PARiOIETERS, I 

I OUTPUT SPECIPICiLTIOU, INITIAL GOTTDITIO?!, BOUNPiRY [ 
; QOU DI TIOUS AND GHiOgfEL JRSAUH SPEOIPIO ATTni\TS__ _ f 


....-"Stage 

-^HYPRO GRAPH IS 
■'•GIVEIT AT U/S 


N 

I 

i GEBIEERATE 

1 PISGHARGE 


1 HYDROGRAPH 


Y ; 





GSUERATE STAGE 
HYPRO GRAPH 


i i 


SET IHITIAE GOHPITIOII 


— 254 


UPPATE TIJffl, RERAIffi YAPUES OF VARIABLES, ' 
ESTIMATE THE APPROXIBJATB Y^iLUES OF Y,Y,Q 
FOR KBIT TIB.® LEVEL 


SETUP ROUTING EQUATIONS 
.INP 

SO LVE T HE MAT.RIX _ 

I 


1 CORREGT THE BSTIMATE P Y , XJM,Q.. 

i- 

^^.^'-"■'CORRECTION^'" 

-■^^--^EERmSSIBLE ERROR-- 





r^' 


kX 

PRINT OUTPUT 


ITERATION > 
§PEOIFIBP VALIpi- 



PRINT ERROR 
JISSSAGB 



FIG. 5.2 FLOW GHA^T OF PROGRM 



27 


In the present study verification has been done by 
the examples worked out by others and comparison of results 
is done among themselves and with them. If two or more 
sets of parameters gives practically identical results, 
one consuming less time is considered the best one. 

It has been observed that the results generally differ 
at the peak. The resultant hydrograph reproducing peak 
values, peak time and shape has been concluded the best. 

In this study the solution is taken as convergent, if it 
converges in 4-0 iteration steps. 



CHAPTER IV 

IlfPLIOII PIHITE ELEM ENT ?!ET H0P 
4 . 1 IFIRODUGIIO H 

In finite element one discretiae the in elements and 
try for solution. The details of the EEM can he seen 
in references^ Cooley and Moin in 1976^^^ 
gave an EEM solution to the flood routing problems. They 
took linear elements and used predictor corrector method, 
and irregular time steps which were less at peaks and more 
else vjhere were used. King in 1977^^^ give another 
solution to the flood routing problem using EEM. He used 
time integration scheme instead of predictor corrector 
method. 

^ ' 2 DEI/JLBE EESGRIPTIOH 

The absic equations were taken in a different form 
here. Eirst and second terms of Eq. 2.1 are gr^'^uped 
together giving new form of continuity equation 

^ - q= 0 (4.1) 

Eq. 2.2 is multiplied by A and Eq. 4.1 by V and by adding 
the results we obtain modified form of momentum equation. 



0(7Q) 
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12 

8 t 


-f 


+ 


OX 


OJ 


’" S^) >- q. — 0 (4.2) 


Equations 4.1 and 4.2 have been used in the finite 
element formulation. 

The system is discretized as used hy Cooley and 
Moin^^^ into linear elements in space for a specified 
time. Hence finite element solution is employed only 
for space variables. King^^^ has suggested the use of 
finite difference time integration scheme for finding the 
values of partial derivatives of variables y and Q with 
respect to time. First describing y as a function of 
time at a particular position in space, 


y=y^ + at + bt°^ ( 4 . 3 ) 

then = a + a b t ( 4 . 4 ) 

Putting the value of t from Eq. 4.4 in Eq. 4.3 we get 


= a + 1 (y - - a a 


0 V 

at t = 0; = a 

then at time /\ 


(®y ) 


~ (yp " ^1) + ^ 


1 


( 4 . 5 ) 
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Here subscript 2 denote the variables at time ^ t. 
Similarly 

^8t ^2 ~ ^.”.”t ^*"^2 (4.6) 

For a = 1 this scheme reduces to the linear integrration 

scheme and for a = 2 the scheme is quadratic integration 
■, . (9) 

scheme. King used the value of a = 1.5. However in the 
present study different values of oc have been tested. 


A formulation based on the method of weighted residuals 
has been used herein to develop the finite element equations. 
The basic Eqs. 4.1 and 4-2 are written in the integral 
form and the shape functions of the finite element approxi~ 
mat ion are used as the weighting parameters. As per 
G-alkerine procedure, 


E r ciS + vlS+oEI +eA3i 


e L. 


2A (S^ - u^ q) dx = 0 


(4-7) 


and /gN 

2 / H. (B ^ + is ^ q) dx =0 (4.8) 


e I, 


5x 


® (e) 

in which is the river reach in the element e and N- 
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are coordinate functions for the node i in the element e. 

The variables Q, y, a, b and are assumed to be 
linearly variable with respect to x in each discrete 
oloment, e.g. 

y = % y^ + y._^^ (4.9) 


In which coordinate fianetion axe given as 


^i = 


and 


F. 


1+1 X 


x: X. 

1 


n+i "■ ^i 


( 4 . 10 ) 


where i is the node bounding the element withing y/hich 

the function is approximated. As variables are assumed 

linear in the element under consider the nodal values 

of elements not bounded by the node will have no effect on 

(o') 

the values of the functions and hence ¥> ^ will be zero 
for all elements not bounded by element e. 

The terms of Sqnsw 4.7 and 4.8 are differentiated 
and integraded as per the indication after putting the 
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different approximating functions. Tlie equations 
for interior nodes will be 

8 Q- ^ ^ n 

«lKl/2) 2 (®l.(1/2) + b+(1/2) 

5t 

c 5 Q- -I 

«l+(1/2) (Q.- Q._^) 

^i+(l/3) " ‘^i^ ^i-(l/3)^^^i’'^i-1^ 

Qi+(1/3) (^i+l ^ ® *'-k-(l/3)^^i^ ^i-1^ 

h+(l/3) ^^i+1 " ^S^^i-(l/2)^i^(l/3)^ou 

^i+(l/2) h+(^/3) ^od ^ 

1 2 
+ ” g (Si^(^/2) h^(^/2) \ ^f i-1 

2 2 

+ (S.^(V2) ^i-(l/4) + ^i+(l/2) V(l/4) ''d ^ 

2 

^fi ^i+(l/2) h+('l/2) ^d ^f i+1 ^ 

’'5(65_^(^/2) "^xn ^i+(l/2) "^xd ^^d^ = ° 

(4.11) 
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and 




( 1 / 2 ) 




■( 1 / 2 ) 


9y 


.Jrrl 

3t 


V(1/2) '-H-(1/4) 


dj. 

at *i+(1/2) ''^l+(1/2) " ‘5i-1^ 

- 6(5i.(i/ 2) + Si^(i/2) 9d> = 0 

( 4 . 12 ) 

where = ± ( 1 ^ 1 + 1 - i^^i) 

^+( 1 / 2 ) 4 - 1 

■^+(1/3) 1 

•%+(l/4) “ ■^+ 1 

I'.P 

Y = Q/A (4-. 13) 


Y, Q, T terms at fractional nodal distance are defined 
in the same v/ay as A has "been defined. The subscript 
u and d indicate the upstream and downstream elements 
with respect to node i. 
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By this , my we get (2 1,^4) equations in 2 H unkno^wns. 
Two more equations are get hy treating the Eq. 4.8 in the 
same way for first and the last elements 


’l + (l/2) ^^1+(i/4) _“_J. + m ''"'2 X 

^ 1/4; — + g-r- ^ 


oy-i 


8 y, 


6 t 


+ 6(Q2 - ^ 5i.(-5/2) <1^) = 0 


(4. -14) 


and 




+ T, 


3 y. 


¥-1 


) 


at ■^^-(V 2 ) " 

+ 6 (Qj^ - Q-n) = ^ ( 4 . 15 ) 


The remaining two equations needed for solution may 
he obtained by the boundary conditions. if the 
upstream stage hydrograph is given then the upstream 
boundary condition will be 

yi ^ fi ('fc^'^^) = 0 (4.16) 

and if discharge hydro graph is given then 


- f^ ( t'^’^'') = 0 (4.17) 

*n » -1 

where ft(t‘ ) carries the same meaning as given in 
Finite difference method. The downstream boundary 
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condition will be 

% ^ = Q ••• ... (4.18) 

where Q is given by Eqs. 2.6 or 2.7 for dynamic and 
Kinematic cases respectively. in the present study 
I'Tornial flow, is calculated by using Manning’s 
equation. 

Kow our equations constitute a set of 2K nonlinear 

algebric equations in 2F unknowns, namely ^^i , ^i , 

3t at 

for i values of 1 to K. However the ^^i and 

\mrwl ■■!««■, iiv 

9t ■ 8t 

can also be furthel? written in the y^ and P®^: 

Eqs. 4.5 and 4.6 respectively. The equations can be 
assembled as follows: 

Rl = (y-i» Qi) =0 

^2 — G''|(y-^5 Q-I ? J2^ ^2^ “ ® 

Rj ~ I'2 (y""!’ ^1’ M2’ *^2’ ^3’ *33^ =0 

^4 — ^2 Qi> ^2’ Q2’ y3’ ^3^ “ 


0 


(4,19) 
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^2i.-1 - 

^4^1-1’ q-1. yi- Qj_. 

^i+1 ’ 

Qi^-|) - 0 

^2i = 

%_i> y^ . Qi 


’^+1 ^ = 0 

>- 

^ V- H- ^ 

>- 



^ 

>— 

- 

■2n-3 = 

%.-2’ 

%T-1 ’ 

Tr %) = 0 

•2n-2 ~ 



yyJ Qjj-) ~ 0 

2bV1 ^ 



= 0 

2n = 

% (yN» 


= 0 


( 4 . 19 ) 


iAlthougli Eq. 4.19 carries 2?T unknowns, yet each 

equation contains a maximum of six unknowns and the 
whole matrix can be put in a banded matrix of width 
seven. The generalized Fewton iteration method is used 
herein to find the solution of this system. Needed 
partial derivatives and the solution may be obtained 
in the same way as described in finite difference 
method. (See Article 3,2). Knowing the values of dy^ 
and dQj_, i having the value from 1 to N, the values of 
yj_ and may be modified^ The method used herein is 
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Gaussian Elimination method using handed matrix technique. 
At the terminal iteration, v;hen values of dy^ and/or 
dQi are less than the permissible errors we take the 
values of variables found as final and advance for time 
stop (j4-2). Eor this first we approximate the values of 
y j_ and equal to (j+1) time level and then find the 
values of time derivatives, in the iterative cycles 
also the values of time derivatives are modified after 
tho completiovi of each iteration. 

4 . 3 PROGRAJ'E'fLITG 

The above PEM method of numerical solution has 
been programed in P0RTRAN>-10 and implemented at DEO- 1090 
system at Indian Institute of Technology, Kanpur. Input 
and output procedure is the same as in Article 3.3. 
Provision for a large number of lateral flows has been 
made. The program can be used for a large number of 
distance and time nodes. Other programming details are 
the same as in Article 3.3. , Plow chart of the program 
is sho?;h in Pig. 3.2. A listing of the PORTRAN" Program 
is included in APPEMDIX - B. 

4.4 GH^ACTERISTIj^S OP PEM^^^^ 

In the present study the selection of distance 
steps has been done by the consideration of stability. 
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and a should be selected such that the program 
is stable, convergent and gives accurate results, in the 
present study it has been found that a = 1.75 satisfies 
the above properties to a large extent. The value of 
cc should be selected between 1.5 2,0 such that it gives 

desired hydrograph. 

In the present study verification has been 
done by the examples worked out by others and comparison 
of results is done among themselves and with them. If 
two or more sets of parameters give identical results, 
one consuming less GHJ time is considered the best one. 
It has been observed that the results generally differ 
at the peak. The resultant hydro graph reproducing peak 
values, peak time and shape has been concluded the best. 
In this study the solution is taken as convergent, if it 
converges in 40 iteration steps. 



OHAPg BR - 7 
_Gqjp>UT'ER STUDIES 

r — jr IT*...*- 

^R^oRMMg e 

Flood 2’outing methods are used to produce the 
relevant outflows by computations. So the accuracy of 
the reproduction may be well defined by a criteria for 
deciding its performance. However, in most of the 
cases coiriplete reproduction is not needed and only certain 
features will be of importance. In the present study as 
the real life floods data were not available an example 
worked out by Yiesmann ^ has been adopted. So 
accuracy will be ijudged from that study only. The peak 
flow, the corre spending stage and the time to peak are 
the most important hydrograph characteristics. However 
the tike to centroid of flow and flow volume may also 
constitute performance criterion. Eventhough effects of certain 
important routing parameters has been studied, but the 
effect of in between computation efforts has not been 
anal.ysed. Further every care has been taken to minimize 
the computational efforts. 

However the performance can not be rated only 
by the accuracy alone but by the convenience of use. 

Gonvenience can be judged by the computer time, simpleness 
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of the model and the range of applicability in the 
practical problems. 

5 • 2 tbsti ug mmLs 

The programs were first tested for programming 

errors. Some attention has been paid to find the limits 

of applicability. 'Phe numerical experimentation is done 

basical]-y for different routing parametrs viz. and 0 

A X 

for finite difference method and and a for finite 

X 

element method. The regular channel geometries were used. 
The following two examples were considered and tested 
with both the methods by using different downstream 
rating conditions. 

5.3 BXAJ.IPLB >- I 

The first example taken is the one as tested by 

Oooley and Moin^^^ on their finite element model. It 

pertains to routing a traiangular discharge hydrograph 

down a rectangular channel and has been adopted from 

f 1 7 "1 

Viessman, et al., \ A rectangular channel 6.1 m wide 
and 3.2 km long carrying a steady uniform discharge of 
23.34 rP/S at 1.83 m depth is subjected to an upstream 
flood wave with a peak of 57 m /S increasing linearly 
in a period of 20 min. This upstream flow decreases 
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-i-inecirly from its peak to 23«34 m^/S in 40 min« 

Muitional properties are = 0.0015 and n = 0,02. 

The explicit method employed by Viessman, et 
utilized a 2 seo. time-step size and a 160 m distance 
step up, Gooley and Moin^^^ used time step of 1 min., 

10 mill, and irreguxar time steps varying between 4 min. 
to 10 min. They kept the dist.-mce interval same. In 
the present study the distance interval taken was 
800 m and time steps tested were 1 min., 2.5 min., 

5 min., and 7.5 min. The problem was solved by both 
finite difference method, and finite element method. 

Both the methods were tested for the Kinematic as 
well as dynamic downstream boundary conditions. The 
resultant hydro graphs have been compared among themselves 
and to those of Kiessman, et al.^*^*^^ and Cooley and Moin^^^. 

In finite difference program the test runs were 
taken for the value of 9 = 0.5, 0.55, 0.60, 0,65 and 
1,00. The program was converging v;ith all the values of 
0 for the time step of 5*00 min. in both the 
kinematic and dynamic cases. However for the value of 
0 = 0.5 - 0.65 the program was also convergent for 
time step of 7,5 min. in the kinematic case. The minimum 
GPU time taken (1.15 seds) was in the kinematic case 
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for a value of 0 = 0,65 » However the results were not so 
^^,ood for time step of 7*5 minutes, it is also a general 
observation tha.t 9 = 0.55 0,60 gives better results 

even though the minimum GFJ time is obtained for 
0 = 0.60 - 0.65. 0 = 1.0 gave much higher CPU 

times. Hence it is thought that 0 = 0,60 may be adopted 
for getting good results as well as low GPU time. Time 
step of 5.0 minutes was adopted for comparison purposes. 

The GHJ time taken is 1,59 sec. for kinematic case and 
1.44 sec. for dynamic case. 

The input and output discharge hydro graphs 
at 1,6 km. and 5.2 km. dovrastream for ^ t = 5 m-in 
and 0 = 0.6 for kinematic and dynamic PPM methods are 

(17) 

plotted in Pig. 5.1. They are compared to the Yiessman^ 

( 

and Gooley and Moin^ ^ hydrographs. The results 
obtained were in between these two methods and kinematic 
and dynamic cases gave similar discharges in this case. ■ 

The input and output stage hydrographs at 5.2 km 
downstream for the above cases are shown in Pig. 5.2. 
Pyuamic models give lower stages and their peak time 
were more when compared to kinematic models. 

In the finite element program the test runs were 
taken for a = 1.0, 1.25, 1.50, 1.75 and 2.0. The GHJ 
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tine for the oonverglng output hydrograph is given in the 
Table 5.1. It is clear from the table that value of 
a does not change the reexcution efforts in the kinematic 



TABLE 5.1 

C.P.U. THE 

FOR PEM 

PROGRM (EXAMPLE l) 


f 

T 



Kinematic 


{ Dynamic 

T 


a 

f A t = 1 . 
I min . 

I 

0 A t = 2.5 
min. 

/a"= ' 

min 

5,0 T/A = 2.' 
. min 

I -r 1 _ _ m r. ■■ 

5 

min. 

1 .00 

2.81 

1.47 

1.05 

2.36 

1.09 

1.25 

2.79 

1.46 

1.05 

2.94 

1.19 

1 .50 

2.74 

1.43 

1.04 

4.10 

1.31 

1.75 

2.70 

1.43 

1.05 

7.20 

1.41 

2.00 

2.67 

1.41 

1.05 


1.56 


case. However in the dynamic case execution efforts 
are more for a higher value of a, GPU time is more in the 
dynamic case. The minimum GPU taken was for ^t = 5.0 min. 
and = 1.50 in the kinematic case. The higher values 
of \y + cl are 0,56, 1.36 and 2,71 for the time 

X 

steps of 1, 2,5 and 5,0 min. respectively and program was 
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stable even for |v i 

I i si values greater than 

The output hvd'”r>xrr'!!iT.v, 

ya-ograph peak shifts to the ri^ht 

with the incroeoino. i bright 

“• The OW time l„ 

the pem model is less th^r. - 

^txat of EDM model. 

Output discharge hydrographs at i.g ^ and 
3.2 kms. downstream for t = 5 min. and c = 1.5 

are shown in Eia- p; 1 ov,^ 

« 5. and are compared to other methods. 

fai elTCs comparable results to other method. 

Kinematic and dynamic models giyes similar results 
tor this problem. 


^ In Pig. 52 stage hydrographs for 1.6 km. downstream 
are given. EEM gave lower stage peak than EDAi. i^namic 
cau. givts lower peak and it occurs with a lag. EEM 
Kinematic model gives good results. 

It is found that if time step is made near zero 
peaks of the discharge hydrographs at 1.6 km. downstream 
computed by EEM for different values of a were found 
to converge at a peak value of 51.4 cumecs. fMs 
value is approximately equal to Yiessman’s^'’'^) results 
for 2 sec. time step explicit solution and hence it 
may be taken as time peak. Eig. 5.5 gives percentage 
results .for higher values of time steps, if they are giv- 
ing stable and convergent results. Dynamic and 
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Kinematic case m^e giving near about same error, 
a - 1 >- 1 ,b gave negative errors; a = 1,75 gave very 

less error; a = 2.0 gave positive errors, a value 
of a = 1.50 - 2.0 may be adopted. In this example 
— 1 ,5 has been taken for comparison purposes while 
a = 2.0 has been taken in next example. Dynamic and 
kinematic cases gave similar errors, 

5,4 gXillFLEy- 11 

In this example hypothetical input hydrograph 
described by the Peason Type III distribution has been 
used. The hydrograph is given by the equation 

Q(t) = Qv, + (Qp - Qi,) (h ) 

where = base flow 

Qp = peak flow 

t„ = time of peak 

P 

P = skewness factor 

O0/7) 

and t« - time of centroid of flow > 

The cross'-sectional area and wetted parameter are given 
by the equations 


1 

^"1 


t:i ■ I; ^ 


= 350 mV sec. 
= 700 mV sec . 

= 300 sec. 

= 1.15 
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A - 0.0008 + 2y^ - 10 y^. + 100 y + 100 

P = 0,0016 y^ + 4y^ ■" 15 J'" +160 

The Manning’s n is 0.025 and bed slope is 
0.0005. The channel is 20 km. long. 

In the present study the distance interval taken 
was 2 kms. Problem was solved by both PDM and PIM. Both 
the methods v/ere tested for the kinematic as well as 
dyniamic down stream boundary conditions. The resultant 
stage and dischc-arge were compared among themselves. 

In the PDM program the test runs were taken fur 
the value of 0 = 0,5, 0,55, 0,60, 0.65 and 1,0. The 
program was tested for ^ t = 21, 25, 40 and 50 mins. 
Smaller values of time could not be tested due to limited 
matrix storage. The minimum computer CFJ time was 2.51 
secs, for /I " 50 min. and 9 = 0.65. The kinematic 
case was not converging for 0 = 1*0 for t > 25 min. 
The dynamic case converged only at t = 21 minutes. 

The GPU time taken is 5.44 secs, for 9 = 0.6, A t > 21 
minutes, dynamic case and 3.17 secs, for 0 = 0.6, 
ilt = 25 min. kinematic case. These resultant 

stage and discharge hydro graphs are plotted in Pigs.5.4 
and 5.5. 
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In the mil Program the test runs were taken 
for a - 1.0, 1 . 25 , 1 . 50 , 1.75 and, 2.0, The GPU time 
foi converging and stable outputs is given in the 
Table 5»2. It can be seen from the table that GPU 

TAJILE 5.2 : GPU Tllffi FOR FEM PROGRAM (EXAMPLE II) 
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min. 
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min* 
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r 

min. 

min. 

min. 

1.00 

25.02 

14.85 

11.07 


26.55 

15.31 

11.64 


1.25 

24.77 

14.59 

10.08 


- 

15.19 

10.53 

- 

1 .50 

24,42 

14.46 

9.52 

11.18 

- 

14.90 

9.86 

11.67 

1.75 

23.69 

14.24 

9.26 

8.80 

■>— 
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9.59 

8.88 
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23 . 2:5 

13.85 

9.10 

7.69 


- 

9.49 

7.22 


time 

is more 

in the dynamic case. The 

minimum CPU 



taken 

was for 

= 

25 min. 

and a = 

= 2.00 

# 



Output stage and discharge hydrographs by FUM 
and FEM at 20 fan. downstream for Jit = 25 min. and 
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and a - 2.0 axe shown in Fig, 5,4 and 5.5 and axe 
compared among themselves, fBM is giving higher stage 
and discharge pea,k8 than FE?J. KinBmatic case gives 
higher stage and discharge peaks. In dynamic case 
stage peak is time lagging. 

It is seen that the program does not converges 
for smaller values of time steps. So proper care should 
he taken in the selection of ^ t and a = 1.5 - 2.0 
may he used for computations. 




diPjaE R YI 

OONGIIJSION 

^ • '• Q OI^OPTJSIOR 

An jj iiiM cilgoicit/luD. feXsnds tlis strong points 

of the two earlier studies has been developed in this 
study* This algonthin has been coinpared. with the 
modified imein's im Scheme. Both the methods were 
tested through numerical experimentation. The following 
signific,ance conclusions are made, 

(l) The proposed PBM which is based on imein's 1970 
work is stable, convergent for 0 = 0,55 - 0.65. 

i 

f) “ 0.6 gives the best results. 

(ii) The proposed REM is found to give stable, convergent 
and accurate results for high time and distance 
steps, FEM program takes less CPU time while 
compared to FBI. It generally have advantage 
over FBM. 

(ill) The value of a = 1.5 - 2.0 gives stable, conver- 
gent and accurate results, A value of a =1.75 
is recommended, 

(iv) The stage hydrographs are more sensitive to the 

downstream boundary conditions while comparing with 



55 


discliargo hydrographs. The peak stage in 
dynamic case is found time lagging. 

^•2 PITT O mP WOIIK 

The proposed FBI has very good advantages like 
use of larger time and distance steps, less computational 
efforts while comparing to other methods. The following 
studies are needed. 

(i) Further studies on. the application to field 
prohlemo, 

(ii) ,P?:U?ametric study regarding selection of a for 
various hydrograph shapes and channel charact- 


eristics. 
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